% 求解三角形(0,0),(1,0),(0,1)上的六个正交基函数系数
% 修复单个变量求解时的访问方式错误
clear; clc; close all;

% 三角形参数
x0 = 1/3;  % 重心x坐标
y0 = 1/3;  % 重心y坐标

% 定义符号变量x和y
syms x y

% 定义三角形区域上的积分函数（符号积分）
% 积分限：x从0到1，y从0到1-x
integrate = @(f) int(int(f, y, 0, 1-x), x, 0, 1);

% 存储基函数系数
num_basis = 6;
coeffs = cell(num_basis, 1);
bases = cell(num_basis, 1);  % 存储基函数表达式

%% 0阶基函数: v0 = a0
% syms a0
% f0 = a0;
% % 归一化条件: ∫∫v0^2 dxdy = 1
% eq0 = integrate(f0^2) == 1;
% % 单个变量求解，返回的是符号变量而非结构体
% sol0 = solve(eq0, a0);
% % 直接对解进行筛选（修复访问方式错误）
% a0_val = double(sol0(sol0 > 0));  % 取正解
coeffs{1} = 1;
bases{1} = 1;

%% 1阶基函数: v1 = a1 + b1*(x-x0)
%syms a1 b1
% f1 = a1 + b1*(x - x0);
% % 正交性条件
% eq1_1 = integrate(bases{1}*f1) == 0;    % 与v0正交
% eq1_2 = integrate(f1^2) == 1;           % 归一化
% % 多个变量求解，返回结构体
% sol1 = solve([eq1_1, eq1_2], [a1, b1]);
coeffs{2} = [0; sqrt(2)];
bases{2} = coeffs{2}(1) + coeffs{2}(2)*(x - x0);

%% 2阶基函数: v2 = a2 + b2*(x-x0) + c2*(y-y0)
syms a2 b2
f2 = a2 + b2*(x - x0)*sqrt(2) + (y - y0)*sqrt(2);
% 正交性条件
eq2_1 = integrate(bases{1}*f2) == 0;    % 与v0正交
eq2_2 = integrate(bases{2}*f2) == 0;    % 与v1正交
%eq2_3 = integrate(f2^2) == 1;           % 归一化
sol2 = solve([eq2_1, eq2_2], [a2, b2]);
coeffs{3} = [double(sol2.a2); double(sol2.b2);sqrt(2)];
bases{3} = coeffs{3}(1) + coeffs{3}(2)*(x - x0)*sqrt(2) + (y - y0)*sqrt(2);

%% 3阶基函数: v3 = a3 + b3*(x-x0) + c3*(y-y0) + d3*(x-x0)^2
syms a3 b3 c3
f3 = a3 + b3*(x - x0)*sqrt(2) + c3*(y - y0)*sqrt(2) + 2*(x - x0)^2;
% 正交性条件
eq3_1 = integrate(bases{1}*f3) == 0;
eq3_2 = integrate(bases{2}*f3) == 0;
eq3_3 = integrate(bases{3}*f3) == 0;
%eq3_4 = integrate(f3^2) == 1;
sol3 = solve([eq3_1, eq3_2, eq3_3], [a3, b3, c3]);
coeffs{4} = [double(sol3.a3); double(sol3.b3); double(sol3.c3);2];
bases{4} = coeffs{4}(1) + coeffs{4}(2)*(x - x0)*sqrt(2) + coeffs{4}(3)*(y - y0)*sqrt(2) + 2*(x - x0)^2;

%% 4阶基函数: 包含交叉项(x-x0)(y-y0)
syms a4 b4 c4 d4
f4 = a4 + b4*(x - x0)*sqrt(2) + c4*(y - y0)*sqrt(2) + d4*2*(x - x0)^2 + 2*(x - x0)*(y - y0);
% 正交性条件
eq4_1 = integrate(bases{1}*f4) == 0;
eq4_2 = integrate(bases{2}*f4) == 0;
eq4_3 = integrate(bases{3}*f4) == 0;
eq4_4 = integrate(bases{4}*f4) == 0;
%eq4_5 = integrate(f4^2) == 1;
sol4 = solve([eq4_1, eq4_2, eq4_3, eq4_4], [a4, b4, c4, d4]);
coeffs{5} = [double(sol4.a4); double(sol4.b4); double(sol4.c4); double(sol4.d4);2];
bases{5} = coeffs{5}(1) + coeffs{5}(2)*(x - x0)*sqrt(2) + coeffs{5}(3)*(y - y0)*sqrt(2) + ...
           coeffs{5}(4)*2*(x - x0)^2 + 2*(x - x0)*(y - y0);

%% 5阶基函数: 包含(y-y0)^2项
syms a5 b5 c5 d5 e5
f5 = a5 + b5*(x - x0)*sqrt(2) + c5*(y - y0)*sqrt(2) + d5*2*(x - x0)^2 + ...
     e5*2*(x - x0)*(y - y0) + 2*(y - y0)^2;
% 正交性条件
eq5_1 = integrate(bases{1}*f5) == 0;
eq5_2 = integrate(bases{2}*f5) == 0;
eq5_3 = integrate(bases{3}*f5) == 0;
eq5_4 = integrate(bases{4}*f5) == 0;
eq5_5 = integrate(bases{5}*f5) == 0;
%eq5_6 = integrate(f5^2) == 1;

vars5 = [a5, b5, c5, d5, e5];
sol5 = solve([eq5_1, eq5_2, eq5_3, eq5_4, eq5_5], vars5);

coeffs{6} = [double(sol5.a5); double(sol5.b5); double(sol5.c5); ...
            double(sol5.d5); double(sol5.e5);2];
bases{6} = coeffs{6}(1) + coeffs{6}(2)*(x - x0)*sqrt(2) + coeffs{6}(3)*(y - y0)*sqrt(2) + ...
           coeffs{6}(4)*2*(x - x0)^2 + coeffs{6}(5)*2*(x - x0)*(y - y0) + ...
           2*(y - y0)^2;

%% 显示结果
disp('三角形(0,0),(1,0),(0,1)上的六个正交基函数系数:');
for i = 1:num_basis
    fprintf('\n%d阶基函数 v%d 的系数:\n', i-1, i-1);
    disp(coeffs{i});
end

%% 显示基函数表达式
fprintf('\n基函数表达式:\n');
fprintf('v0 = %.4f\n', bases{1});
fprintf('v1 = %.4f + %.4f*(x-1/3)\n', coeffs{2}(1), coeffs{2}(2));
fprintf('v2 = %.4f + %.4f*(x-1/3) + %.4f*(y-1/3)\n', coeffs{3}(1), coeffs{3}(2), coeffs{3}(3));
fprintf('v3 = %.4f + %.4f*(x-1/3) + %.4f*(y-1/3) + %.4f*(x-1/3)^2\n', ...
        coeffs{4}(1), coeffs{4}(2), coeffs{4}(3), coeffs{4}(4));
fprintf('v4 = %.4f + %.4f*(x-1/3) + %.4f*(y-1/3) + %.4f*(x-1/3)^2 + %.4f*(x-1/3)(y-1/3)\n', ...
        coeffs{5}(1), coeffs{5}(2), coeffs{5}(3), coeffs{5}(4), coeffs{5}(5));
fprintf('v5 = %.4f + %.4f*(x-1/3) + %.4f*(y-1/3) + %.4f*(x-1/3)^2 + %.4f*(x-1/3)(y-1/3) + %.4f*(y-1/3)^2\n', ...
        coeffs{6}(1), coeffs{6}(2), coeffs{6}(3), coeffs{6}(4), coeffs{6}(5), coeffs{6}(6));

%% 验证正交性
fprintf('\n正交性验证（应接近0或1）:\n');
for i = 1:num_basis
    for j = 1:i
        ortho = double(integrate(bases{i}*bases{j}));
        fprintf('∫v%d*v%d dxdy = %.6f  ', i-1, j-1, ortho);
    end
    fprintf('\n');
end

%% 高斯积分函数 - 在三角形上积分
function val = gauss_integrate(f, vertices,x,y)
    % 在三角形上使用高斯积分计算函数f的积分
    % 三角形高斯积分点(3点公式)
    %gauss_points = [1/6, 1/6; 2/3, 1/6; 1/6, 2/3];
    %weights = [1/6, 1/6, 1/6];  % 权重
    gauss_points = [0.1012865073235,0.1012865073235;
                    0.7974269853531,0.1012865073235;
                    0.1012865073235,0.7974269853531;
                    0.4701420641051,0.0597158717898;
                    0.4701420641051,0.4701420641051;
                    0.0597158717898,0.4701420641051;
                    0.3333333333333,0.3333333333333];
    weights = [0.1259391805448;0.1259391805448;0.1259391805448;
        0.1323941527885;0.1323941527885;0.1323941527885;0.225];
    % 计算三角形面积
    area = 0.5 * abs(det([vertices(2,:)-vertices(1,:); vertices(3,:)-vertices(1,:)]));
    
    val = 0;
    %syms x y
    f_handle = matlabFunction(f, 'Vars', {x, y});
    
    % 映射到实际三角形
    for i = 1:3
        xi = gauss_points(i,1);
        eta = gauss_points(i,2);
        
        % 从参考三角形到实际三角形的映射
        x_coord = vertices(1,1) + xi*(vertices(2,1)-vertices(1,1)) + eta*(vertices(3,1)-vertices(1,1));
        y_coord = vertices(1,2) + xi*(vertices(2,2)-vertices(1,2)) + eta*(vertices(3,2)-vertices(1,2));
        
        % 计算函数值并累加
        val = val + weights(i) * f_handle(x_coord, y_coord);
    end
    
    % 乘以面积得到最终积分值
    val = val * area;
end